Computing Dynamic Output Feedback Laws 

Jan Verschelde* Yusong Wang"'" 



2 June 2003 



Abstract 

The pole placement problem asks to find laws to feed the output of a 
plant governed by a linear system of differential equations back to the 
input of the plant so that the resulting closed-loop system has a desired 
set of eigenvalues. Converting this problem into a question of enumerative 
geometry, efficient numerical homotopy algorithms to solve this problem 
for general Multi-Input-Multi-Output (MIMO) systems have been pro- 
posed recently. While dynamic feedback laws offer a wider range of use, 
the realization of the output of the numerical homotopies as a machine to 
control the plant in the time domain has not been addressed before. In 
this paper we present symbolic-numeric algorithms to turn the solution 
to the question of enumerative geometry into a useful control feedback 
machine. We report on numerical experiments with our publicly available 
software and illustrate its application on various control problems from 
the literature. 
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1 Introduction 

Given a linear system of differential equations and a list of eigenvalues, the pole 
placement problem asks to find laws to feed the output back to the input so 
that the resulting closed-loop system has the same eigenvalues as the given list. 
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While executing a static feedback law corresponds to a simple matrix-vector 
multiplication, a dynamic compensator has several internal states. 

The theoretical solution of this problem was found in the Schubert calculus, 
see jjU for the static, and |251 1^| for dynamic compensators. Because of 
its importance to practical applications, the development of algorithms for this 
problem was stated as an open problem |2H] (see also jS])- The first homotopy 
algorithms were proposed in |14|. defining a numerical Schubert calculus. The 
Pieri homotopy algorithms of |14| were improved and generalized to dynamic 
feedback in ^5]. In [201, the numerical performance of these homotopies for 
static feedback was improved. 

Our current paper is a sequel to our previous paper as a continuing effort 
to apply the homotopy algorithms to practical applications of linear systems 
control. Rather than generating random matrices as inputs to our algorithms 
we look to build case studies with examples from the literature. In we 
applied static output feedback to use pole placement to keep a satellite in orbit. 
In this current paper we consider the application of dynamic feedback laws to 
this and other applications. See ^B] for the relation of the pole placement 
problem with inverse eigenvalue and matrix extension problem. 

As expected, with the more general dynamic feedback laws we can cover 
a wider range of applications, but at the same time we experienced the real 
need for more advanced algorithms from numerical algebraic geometry to deal 
better with over and underdetermined problems. We define these mathematical 
problems in the next section. 

The heart of this paper concerns the realization of dynamic compensators. 
As these compensators are computed in the frequency domain, at the end of 
the paths defined by the homotopies, any application demands the description 
of the controllers in the time domain as a tuple of matrices. While realization 
algorithms are covered in standard graduate textbooks (e.g., pQ and our 
output is numerical, i.e.: subject to roundoff and known only with limited ac- 
curacy. Reports on numerical implementations of algorithms which manipulate 
polynomial matrices are scarce, we found only one paper |12| . which presents 
a numerical study on a commercial implementation, but leaves the implemen- 
tation of a numerical Smith normal form as an open question. Note that the 
Smith normal form plays an important role to obtain a minimal realization of 
the transfer function of Multi-Input-Multi-Output (MIMO) systems. 

Section three starts with the description of an algorithm to compute the 
greatest common divisor of two polynomials with approximate coefficients, as 
this is the cornerstone function in the implementation of a numerical Smith 
normal form. As we experienced and illustrate, the naive application of the 
Euclidean algorithm for an approximate GCD can fail miserably. That many 
algorithms to solve problems with exact data turn out to be numerically unstable 
is a growing concern in computer algebra, which has led to hybrid symbolic- 
numeric computation |SJ. In particular, the approximate GCD problem has 
received a lot of research attention in recent years, see for instance [2], |E|, EH, 
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P21» El, EH and EB1 

Our approach is to consider the input polynomials not as polynomials with 
approximate coefficients, but as polynomials defined by approximate roots. So 
we view polynomials not symbolically (as sums of monomials with approxi- 
mate coefficients) but geometrically (as the collection of roots in the com- 
plex plane, known with limited accuracy). The idea to computing the GCD 
by matching common approximate roots (within a certain tolerance) can be 
found in Compared (see [21]) to alternative constructions of approxi- 

mate GCDs, this method works well for input polynomial with multiple or 
clustered zeros. Algorithmically, we can immediately apply any numerical root 
finding method (see |23j for a recent survey). To compute the polynomials 
in the extended GCD problem (i.e.: the k{x) and l(x) in the combination 
k{x)a(x) + l(x)b(x) — GCD(a(x), b(x))), we apply Newton interpolation. 

As usual, the link between algorithms and applications is software. We 
describe an extension to the software package PHCpack [[2] in section five before 
illustrating its application on various control problems from the literature. 

2 Problem Statement 

The algorithmic framework is defined by the transitions between the time and 
the frequency domain (which leads to the symbolic- numeric calculations) . Math- 
ematically, we count the number of conditions to determine the dimension of 
our solution space. 

In the time domain, we picture the control of a plant with a dynamic feedback 
law in Figure ^ 
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Figure 1: Control of an m-input and p-output plant by a qXh order dynamic 
compensator in the time domain. 



The machine we want to control is given by three matrices (A, B, C), defining 
a system of first-order differential equations in the time domain. The input of 
the Pieri homotopies is sampled in the frequency domain. So to apply these 
homotopies, we need to translate the input from the time to frequency domain 
and to realize the output as a tuple of matrices. Tabled defines the three stages 
in the data flow. The labels on the arrows in Figure[3correspond to the numbers 
in Tabled 
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(1) 


With plain linear algebra we compute the input for the Pieri 
homotopies, sampling points from the plant we wish to control. 


(2) 


Pieri homotopies compute solution maps of degree q, 
which are then brought into another format. 


(3) 


Given the transfer function of the compensator, we realize 
the compensator by a tuple of four matrices (F, G, H, K) . 



Tabic 1: Transitions in applying the Pieri homotopies. 



Our problem is thus to process the output of the Pieri homotopy algorithms 
and to apply realization algorithms to the numerical output. 
In |15|. the following equivalence relation was derived: 
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Equation (JTJ is the characteristic equation of the closed-loop system. Via el- 
ementary row and column operations, this equation can be rewritten into (0, 
which separates the given data (A, B, C) from the unknown (F, G, H, K). Equa- 
tion @ shows the geometric problem: we are looking for curves which produce 
p-planes in C m+P which meet given m-planes sampled at prescribed values for s. 
In Figure|21we show the transition between the time and the frequency domain. 



x = Ax + Bu 
y = Cx 



z = Fz + Gy 
u = Hz + Ky 



C(sI n -A)- l B 



(2) 



(3) 



Hislq-F^G + K 



Figure 2: Transitions between time and frequency domain, as defined in Tabled 

When controlling a machine with n internal states with a controller using q 
internal states, we can place n + q poles. On the other hand, the dimension of 
the geometric problem is mp + q(m + p). Depending on values for n, m, p, and 
q, we distinguish three cases: 

1) n + q < mp + q(m + p) underdetermined: For a generic machine, there is 
a set of feedback laws. The set has dimension mp + q(m + p) — n — q, and 
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has degree d(m,p,q), i.e: for a generic choice of the parameters, we have 
d(m,p,q) complex feedback laws. This degree can easily be computed 
combinatorially. 

2) n + q = mp + q(m + p) dimension zero: For a generic machine, there are 

exactly d(m,p,q) complex feedback laws. Every feedback law places all 
n + q poles at the desired locations. It may be that no feedback law has 
all its coefficients real, see pi IIP). 

3) n + q > mp + q(m + p) overdetermined: For a generic machine, there are 

no feedback laws which place all n + q poles at the desired locations. 

The numbers m, p, n are fixed, given on input. We can choose q to arrive 
always in a favorable condition and get feedback laws. In particular, if a system 
is overdetermined for static output feedback (n > mp), we can always find q 
to make it convert to underdetermined or dimension zero case with dynamic 
feedback. This is because q(m + p) always increases faster than q. 

Presently, we resolve the undetermined case by choosing additional input 
planes to the geometric problem. For nongeneric machines, the dimension of the 
set of feedback laws is higher than expected. Recent advances with homotopics 
(see e.g. 23 ESI EH EH) allow to treat positive dimensional solution sets, but 
we defer the application of these recent homotopy methods to a future paper. 



3 Symbolic-Numeric Calculations 



To transform a solution obtained by the homotopy methods into a description of 
a machine to execute the dynamic feedback law, we need to calculate the Smith 
normal form to compute the inverse of a matrix with polynomial entries. More 
precisely, the output of the homotopies as an (m+p)-hy-p matrix of polynomials 
in s: 

- U(s) 
V(s) 



(3) 



where U{s) is a p-by-p matrix and V(s) is an m-by-p matrix of polynomials in 
s, satisfies 

U(s) C(sl n - A)^ 1 B 

V(S) l m 

for the given poles. We can right multiply (0J by 
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The result of this multiplication is 
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Since the multiplier matrix JSJ is of full rank, its determinant is nonzero and the 
original intersection condition remains. This multiplication does not affect the 
input conditions, which are at the right part of @- By comparing J2J with ©, 
we can apply the realization algorithms to extract (F, G, H, K) from the matrix 

y( s )c/- 1 (s). 

In the next subsection we show how the calculation of a Smith normal form 
requires the calculation of greatest common divisor. 



3.1 Numerical Smith Normal Form 



The Smith normal form of a matrix of polynomials can be seen as the discrete 
analogue of the singular value decomposition. More precisely, for any n-by- 
m matrix A(x) whose entries are polynomials in x, there exists a unimodular 
n-by-n matrix P(x) and a unimodular m-by-m matrix Q(x) so that 

P(x)A(x)Q(x) = D(x), (7) 

where D(x) is an n-hy-m matrix which has only nonzero polynomials in x on 
its diagonal. Furthermore, denoting the ith element on the diagonal of D(x) by 
D il we have that Di divides D i+ i. The matrix D(x) is called the Smith normal 
form of A(x). Since unimodular matrices are invertible, we can rewrite A(x) as 

A(x) = P- 1 (x)D(x)Q- 1 (x), (8) 

which reveals the following expression for the inverse of A(x): 

A~\x) = Q(x)D- 1 (x)P(x), (9) 

which of course only exists if D(x) has full rank. We get the inverse D~ 1 (x) by 
inverting every entry on the diagonal of D{x). 

The Smith normal form can be computed by solving the extended greatest 
common divisor problem. In particular, we wish to find polynomials k(x) and 
l(x) so that the following equality holds: 

d(x) = GCD(a(x),b(x)) = k(x)a(x) + l(x)b(x). (10) 

The calculation of the greatest common divisor is used to reduce columns 
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(12) 



The matrices we use in the column and row reductions have determinant equal 
to one, they are called unimodular. 

Collecting the column reductions in P(x) and the row reductions in Q(a;), 
we can reduce any polynomial matrix A(x) to a diagonal form D(x), and so 
obtain the Smith normal form of A(x). 
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3.2 Numerical Greatest Common Divisor 



Suppose we are given two polynomials in one variable and with complex coef- 
ficients, only known approximately. We denote the given polynomials by a(x) 
and b(x), and d(x) = GCD(a(x), 6(x)), with deg(d(x)) = r. 

If we apply the algorithm taught in elementary school to compute the great- 
est common divisor of two natural numbers, then we repeatedly have to divide 
polynomials. This repetitive division is numerically unstable as it involves the 
subtraction of polynomials, which may lead to a dramatic loss of significant dig- 
its when the polynomials have coefficients of equal magnitude. Our approach 
to computing the greatest common divisor of polynomials with approximate 
coefficients hinges on two operations: root finding and interpolation. Any in- 
troductory course in numerical analysis describes numerically stable algorithms 
for these two operations. 

For a given tolerance e > 0, we define the numerical greatest common di- 
visor of a(x) and b(x) as the monic polynomial whose roots are common to 
a(x) and b(x) within the given tolerance e. More precisely, if a(a^) = 0, for 
i = I, 2, . . . , deg(a(x)) and = 0, for i — 1, 2, . . . , deg(6(x)), then we can 

rearrange the indices of the roots of the two polynomials so that the r common 
roots appear with lowest indices. Then we write 



r deg(o(x)) dcg(a(>)) 

a ( x ) = T\.( x ~ a ^ IT {x - a{) = di{x) Y\ (x-tti) (13) 

i—l i—r+1 i— r+1 

and 

r deg(6(a:)) deg(6(x)) 

6(x) = JJ(x-A) II (x - I3i) = d 2 (x) J] (x-ft) (14) 

i—l i— r+1 i— r+1 

where \on— /3j| < e, for i = 1, 2, . . . , r, and \oti — (3j\ > e, for all i and j with index 
higher than r. The polynomials di(x) and e?2(x) are numerical approximations 
for the greatest common divisor d(x) of a(x) and b(x). 

Now that we have numerical approximations for d{x), we want to find k(x) 
and l(x) defined in (|1U|I . We determine k{x) by interpolation at those roots of 
b{x) not shared by a(x) replacing in (fLU|) x by for i = r + 1, . . . , deg(6(x)): 

d(A) = fc(AMA) or fc(A) = ^T, fori>r. (15) 

Note that as i > r: a((3i) ^ 0. The interpolation conditions in l|15|) determine 
k(x) uniquely as a polynomial of degree deg(6(x)) — r — 1. Analogously, we de- 
termine ^(x) by interpolation at those roots of a(x) not shared by &(x) replacing 
in the equality (fTTTfl i by a.;, for i = r + 1, . . . , deg(a(x)): 

d(ai) = l(ai)b(ai) or Z(a;) = ^^4, for i > r. (16) 

b(ai) 

Note that as i > r: b(cti) ^ 0. The interpolation conditions in 1|16|) determine 
^(x) uniquely as a polynomial of degree deg(a(x)) — r — 1. 
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3.3 Numerical Experiments 



The algorithms described above have been implemented in C. In this section we 
list some numerical results, obtained by practical comparisons between our new 
algorithm and the elementary approach, for random and specific input data. 

We call the elementary school algorithm by repetitive division the "naive al- 
gorithm" and the algorithm with root finding the "advanced algorithm" . When 
the tolerance of 1CP 8 is used to decide whether two numbers are equal, the naive 
algorithm runs much faster than the other one, while the advanced algorithm 
is more numerically stable when the polynomial degree is less than 30 and the 
degree of the GCD is less than 15. Although the advanced algorithm is rela- 
tively slow, the time spent for each GCD computation is just trivial, less than 10 
milliseconds even for the polynomial of degree 30 and the GCD of degree 15. 
Table [3 summarizes this experiment. 



degree of a and b 


degree of GCD 


naive (%) 


advanced(%) 


5 


3 


100.0 


100.0 


10 


5 


99.2 


100.0 


15 


8 


99.5 


100.0 


20 


10 


99.5 


100.0 


25 


13 


98.2 


100.0 


30 


15 


88.2 


100.0 



Table 2: Comparison of the naive algorithm and the advanced algorithm on 
random data. The first column lists the degree of the given polynomials a and b. 
The second column lists the degree of the GCD. The data in columns three and 
four shows the percentage of success for 1000 tests of the two algorithms. 



We also compared the performance of the algorithms for some specific data. 

In the first experiment, the leading coefficient of the polynomial b is 10 -5 
times leading coefficient of the polynomial a. We found that the naive algorithm 
fails completely when the degree e of the input polynomials and the degree r of 
the GCD satisfy the following relation: 

_ J 2r — 1 for r is odd , 
1 2r for r is even 

Other relations between the degrees work fine. While the advanced algorithm 
remains numerically stable no matter what the relations between the degrees are. 
See the column with header specific numbers(l) in Table [3| for the experiment 
of this case. 

Second, we test with some other specific numbers, say, if the higher degree 
coefficients of two polynomials are very near each other, then the two algorithms 
perform totally different. In particular, we tested polynomials with equal lead- 
ing coefficients and with difference between the second-highest coefficient equal 
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to 1CP 5 . When the degree of the input polynomials is smaller than 5 and the 
degree of the GCD is smaller than 3, both algorithms work perfectly. When 
the degree of the input polynomials is larger than 10 and the degree of GCD is 
larger than 5, the naive algorithm fails completely, but the advanced algorithm 
shows the same numerical stability as on the random number case. See columns 
with header specific numbers(2) in Table |3| for the experiment of this case. 





specific numbers(l) 


specific numbers(2) 


degree 


degree 


naive 


advanced 


naive 


advanced 


of a and b 


of GCD 


method(%) 


method(%) 


method(%) 


method(%) 


5 


3 


0.0 


100.0 


99.9 


100.0 


10 


5 


99.1 


100.0 


0.1 


100.0 


15 


8 


98.8 


100.0 


0.0 


100.0 


20 


10 


0.3 


100.0 


0.0 


100.0 



Table 3: Comparison of the naive algorithm and the advanced algorithm on 
specific data. The first column lists the degree of the given polynomials a and b. 
The second column contains the degree of the GCD. The data in columns three 
to six shows the percentage of success for 1000 tests of the two algorithms. 



These experiments provide us with practical evidence that for both random 
and specific inputs, the advanced algorithm shows its strong numerical stability. 
Concerning the speed of the GCD algorithm, the time needed of the realization 
algorithm is negligible compared to the calculation time of the feedback laws 
with homotopies. 

4 Realization of Multi-Input Multi-Output Sys- 
tems 

In the previous section, we gave the derivation of the transfer function of the 
dynamic compensator as the output of the homotopies. We will give a modified 
algorithm based on ^ pages 389-416] to obtain minimal or irreducible realiza- 
tions, which realize a system with the least number of dynamic elements. These 
modifications were made to fit the output format of the software used to com- 
pute the feedback laws. The necessity of the modifications will be discussed at 
the end of this section. We will show how to obtain realizations {F c , G c , H c , K c } 
of the transfer function T(s) in controller form first. Then, we will use the prop- 
erty of the output of homotopies to show the realizations are irreducible, so they 
are also controllable. 

From JSJl, the transfer function can be written as 

T(s) = V(s)U-\s). (18) 
In accordance with convention, wc would replace V(s) by N(s) which stands 
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for numerator, replace U(s) by D(s) which stands for denominator. 

According to Theorem 3.3 in ^ page 391], realizations exist if and only if 
T(s) is a matrix of rational functions and satisfies 



lim T(s) < oo, 



(19) 



i.e., if and only if T(s) is a proper rational matrix. Given the transfer function 
matrix T(s) = N(s)D~ 1 (s) as a (m x p) proper rational matrix. Let dj = the 
highest degree of jth column in the D(s) (dj > 0,j = 1, 2, . . . ,p). Define 



A( s )=diag(s dl ,..., S ^), 



and 



S(s) = block diag 



/ 


1 

s 


V 





3 = 1, 



(20) 



(21) 



/ 



If dj — 0, just skip that column and continue to fill the next column of the S(s) 
matrix. Note that S(s) is an q(= X)j=i dj) xp polynomial matrix. Write 

D(s) = D h A(s) + DiS(s) (22) 

Dh is the highest column degree coefficient matrix of D(s). For example, if 
3s 2 + 1 2s 1 

, then the highest column degree coefficient matrix 



D(s 



2s 



3 2 
1 



, and DiS(s) given in (|22() accounts for the remaining lower 

column degree terms D(s), with Di being a matrix of coefficients. 
In general, \Dh\ ^ 0, and define p x p and p x q matrices 



respectively. Then F Cl G c can be determined from 

F c = F c + G c F p , G c — G C G P , 

where F c = block diag^i, F% . . . , F p ] with 









G c = block diag 












/ 


' " 







V 


1 







e R 



eR d >,j = l,. 



. ,p 



(23) 
(24) 

(25) 
(26) 
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When dj = 0, we just skip the corresponding Fj and continue to fill the F c 
matrix with the Fj+i matrix; we also need to add a zero column at the jth 
column of the G c matrix. 

Then we can determine H c and K c such that 

N(s) = H c S(s) + K c D(s), (27) 

and note that 

K c = lim T(s). (28) 

s — >oo 

Therefore, H c can be determined from (|27|l . 

An qth-order realization of T(s) in controller form is now given by the equa- 
tions 

z c = F c z c + G c y, u = H c z c + K c y. (29) 

According to the format of the output of the software, q = Y^j=i dj is equal 
to the minimal order of the dynamic compensator. Therefore, this algorithm 
gives us a minimal realization of the transfer function matrix T(s) and the result 
is also observable. 

The main difference between the modified algorithm and the original algo- 
rithm given in ^ is that the original algorithm limits dj > 1, while our modified 
algorithm works for dj > 0, where dj is the highest column degree of jth col- 
umn in the D(s). Some dj must be equal to zero when the number of output is 
larger than the order of the dynamic compensator. In this case the modified al- 
gorithm become necessary. The correctness of the modified algorithm is verified 
with experiments. 



5 Software 

The dynamic feedback laws were calculated with the aid of PHCpack p2] . While 
the second public release of PHCpack implemented the dynamic pole placement 
problem in its geometric form, additional software had to be written, concerning: 

0. a limit on the number of feedback laws; 

1. an interface between Ada and C; and 

2. a collection of C routines for the realization. 

The limit on the number of feedback laws was imposed as a matter of conve- 
nience, to control the practical complexity. We elaborate the other two items 
in the following subsections. 
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5.1 AC interface to PHCpack 



PHCpack is written in Ada, while the programs to process the feedback laws 
are in the lower level language C. 

We can build a portable interface to the Ada routines in PHCpack with 
C functions because the language Ada has the pragma Import construction to 
call routines from other languages such as C and it supports conversions for 
C integers, doubles, and arrays of these C types. Furthermore, the gnu-ada 
compiler provides a mechanism to call Ada routines from a C main program 
and to call C functions from Ada. As the gnu-ada compiler is integrated in 
the gcc compilation system, our interface is portable. In particular, we ran our 
implementation successfully on SUN workstations running Solaris and on PCs 
running Linux and Windows. 

To exchange data efficiently, programs in Ada or C should define exchange 
protocols of structured data types into basic data types for which automatic 
conversions are supported. More precisely, we represent structured data types 
into arrays of doubles and arrays of integers. The language C is restricted 
in returning dynamically allocated variables. Therefore, data allocated in a 
C function is passed by the C function calling an Ada function for further 
processing of the data. 

A typical sequence of calls goes as follows. First a C function gathers problem 
data and prepares the input to an Ada routine of PHCpack. The Ada routine, 
called from C, uses path tracking to solve the problem, and then calls a C func- 
tion to process the results obtained with PHCpack. So the C programmer who 
uses PHCpack should thus provide two C functions: one to prepare the input 
and one to process the output. This "hand-in-glove" interface is appropriate for 
a C programmer collaborating with an Ada programmer (which is the case of 
the authors), who only have to agree on the prototypes of the routines. 

5.2 The Organization of the Software 

In Figure |3 the arrows indicate the order of function calls in the computation 
and realization of the dynamic feedback laws. In this section we give a short 
description for each of the procedures, some written in C, others in Ada. 

C to Ada Ada to C 

ts_feedback(C) j =>^pieri_solver(Ada) ^) 5^ pieri_sols(C)J 




( realization(C) ~) a»(test_results(C) ) a-Tverify(MATLAB2) 



Figure 3: Organization of the Software 
ts_feedback(C): ts_feedback reads all the input information from a file, includ- 
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ing the number of the internal states n, the input dimension m, the output 
dimension p, the number of the internal states for the dynamic compen- 
sator q and the number of output feedback. Also the user should give the 
A, B. C matrices (or let the matrices be generated randomly) of the given 
plant and n + q eigenvalues. tsjeedback.c computes C(sl n — A)~ l B at 
the interpolation points as the input planes. With C to Ada interface, 
we pass the arrays of the input planes and the interpolation points to the 
pierLsolver (an Ada procedure in PHCpack). 

pierLsolver ( Ada) : The pierLsolver calculates the corresponding dynamic out- 
put feedbacks and passes them to the C program pieri_sols.c. 

pierLsols(C): With Ada to C interface, the arrays in Ada form are converted 
to the form in C. Then pierLsols.c calls realization. c and tests the results. 

realization(C): We use modified realization algorithm based on to get the 
realization of the dynamic output feedback, organized as follows: 

(a) Get the transfer function T(s) = N(s)D(s)^ 1 from the output of 
the Ada program. The inverse of a polynomial matrix is a rational 
polynomial matrix and it is mainly done by Poly_Smith. See (10) 
in section 3 for more detail about how Smith form can work for the 
inverse of a polynomial matrix. 

(b) The realization function implements the modified realization algo- 
rithm to get a minimal realization of the dynamic compensator. 

(c) Evaluate the transfer function T(s) at some random point and com- 
pare it with the result after realization (H(sl q — F)~ 1 G + K) at the 
same point. If the values are the same, the realization is correct. 

test_results(C): We can evaluate equation (j2J at the given poles and calcu- 
late the determinant with the previous result. If the determinant is zero, 
the pole is the eigenvalue of the closed-loop system. As mentioned above, 
equation J5J) is algebraically equal to equation (JJJ , which is the character- 
istic equation of the closed-loop system. 

verify (MATLAB): Finally, a MATLAB script verifies the results by compar- 
ing the computed poles with the given poles and finding the condition 
number for each given pole. 



5.3 Availability of the Software 

The C routines for the numerical realization algorithms written by Yusong 
Wang are available at http://www.math.uic.edu/~jan in the distribution of 
release 2.2 of the source code of PHCpack (see (31], the first version of this 
package is archived by Netlib). In particular, the collection of routines can be 
found in the directory "Feedback" of the source code. Also the input data files 
for the applications discussed in the next section are available in this directory. 
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6 Applications 



In this section we illustrate the usefulness of our approach. Some application 
examples will show it is necessary to find some dynamic output feedback laws 
for some specific situations, especially for the cases there is no feedback laws or 
no real feedback laws at some given poles for the static feedback problem. 

We assume that the input data for our applications is given as (A,B,C), 
i.e.: a triple of three matrices of the linear system x = Ax + Bu, and y = Cx, 
where x, u, and y are vectors of internal states, input, and output respectively. 



6.1 Satellite Trajectory Control 



This application concerns the design of output feedback laws to keep a satellite 
in orbit. We treated this problem in our previous paper with static output 
feedback and with ad hoc MATLAB scripts to aid the manipulations. In |35|. 
we already found two real output feedback laws for static case. Now we want 
to get some feedback laws for dynamic case as an extension of our method. 

The linearized state-space equations for the satellite problem are given below: 



A = 



B 





0.3578 





1.3411 





1 



-0.5596 





1.0867 





0.8525 
1 




(30) 



We choose value for each parameter randomly. While this may seem unrealistic, 
our choice can be justified by an appropriate selection of units. More detail 
about the model can be found in [7| I16| . We define C as some random matrix 
- which can be interpreted as a random projection of the states onto a plane. 

For this satellite example, we have n — 4, m = p = 2, which is dimension 
zero for the static case. When q = 1 , the system become underdetermined since 
n+q = 5 < mp+q(m+p) = 8 and there are three degrees of freedom. We choose 
the eigenvalues as {=^, -5, -7, -3.0, -0.1068, -0.7834, -0.9582). The 

last three of eigenvalues are randomly selected. We found two real feedback 
laws and six complex complex feedback laws. We also use the MATLAB script 
to verify the results and calculate the condition number for the eigenvalue by 

1 . . 

— jj r, %= 1,2,. ..,71+0, (31) 

|y; r x *l 

where vectors Xj and y^ denote the unit right and left eigenvectors of the closed- 
loop system. See ^] page 323] for the derivation of |yfxi| as the reciprocal of 
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the condition number for the eigenvalue A^. By substituting the result into the 
closed-loop system, we find the relative difference of the computed eigenvalues 
vs. the given eigenvalues is bounded by 10 -11 and the order of condition num- 
bers are at most 10 3 , calculated with 1)31(1. The total CPU time spent on the 
calculation of the dynamic feedback laws is 2 seconds and 320 milliseconds on 
a 2.4GHz workstation running Linux System. 



6.2 Numerical Example 

In this section we report on two numerical examples in the literature |27| |36| . 

Numeric Example A: We will use the example in to illustrate the follow- 
ing situation: when a system is overdetermined for static output feedback 
(n > nip, q = 0), for which no feedback laws can be found at the desired 
poles, we can choose a q to make it convert to underdetcrmincd or dimen- 



sion zero case (n 
laws. 



q < mp + q(m + p)), for which we can find feedback 



Consider the system, given by 



A = 





1 
1 













1 





0-101 



B 



c 



1 
1 



(32) 



Here, m = p = 2, n = 6, therefore n > mp which is overdetermined for 
static case. We will choose a q to make the system underdetermined. The 
minimum possible q is 1, then we have n + q = 7 < mp + q(m + p) = 
8, so there is one degree of freedom. With the calculation of PHCpack 
and the C interface, we can easily get some dynamic compensators of 
McMillan degree q — 1 to control the system. We verified the result by 
comparing the given and computed eigenvalues with a MATLAB script. 
For a choice (-0.1, -1.5, -0.9, -0.7, -6.0, -3.5, -8.0) of 7 eigenvalues and 
an additional pole —0.1053 which is generated by random, the relative 
difference of the computed eigenvalues vs. the given eigenvalues of the 
closed- loop system is bounded by 10" 9 . In this specific example, we find 
8 solutions and 4 of them are real (sometimes 6 are real, depending on the 
additional input plane which are randomly generated) . The order of the 
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condition number computed with (|31[l is no more than 10 3 for all of the 
given eigenvalues. The total CPU time spent for this numeric example is 
around 3 seconds with the workstation mentioned above. 

Numeric Example B: The second example can be found in 36, Example 3.7], 
with n = 8 and m = p = 3. 

The system is given as below 



.4 






-1 

















1 


1 


2 








1 








-2 





-1 








5 

















1 





-7 








-2 





-1 





1 


4 








2 





-2 








2 








3 














-1 


1 





-2 





-1 








1 





1 


-1 



D 



2 
1 

-3 
1 
4 
5 
1 

-1 



(33) 



C = 



10 





-2 
10 
1 



This system is underdetermined for static output feedback (n < mp, q = 
0). When the given poles are (-0.8090+0.5878^, -0.9511-0.3090i, -0.3090- 
0.951H, -0.3090+0.9511i, -0.9511+0.3090i, -0.8090-0. 5878i, -0.5878+ 
0.8090i, -0.5878 - 0.8090i, -0.1883), in which the first 8 are picked at the 
unit circle and the last one is a negative real number chosen by random, 
we find 42 feedback laws with our software and 4 of them are real. The 
relative difference of the computed eigenvalues vs. the given eigenvalues is 
10~ 9 . The order of the condition number is bounded by 10 5 if computed 
with 131|) . The total CPU time spent on this example is around 50 seconds 
to find all 42 feedback laws. It only takes 1 second and 180 milliseconds 
if user just needs one feedback law. 

6.3 Aircraft Control 



It may be that for a given selection of poles, all static feedback laws have 
coefficients with nonzero imaginary parts. In this case, we will design dynamic 
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output feedback laws, and exploit the additional freedom to place the poles at 
their originally selected locations. The model from MathWorks [21] could be 
used to illustrate this kind of situation. It also can be found in any standard 
text in aviation for a more complete discussion of the physics behind aircraft 
flight. The state-space equations are given below. 

The jet model during cruise flight at MACH=0.8 and H=40,000ft. is 



A = 



B 



-0.0558 
0.5980 
-3.0500 


0.0073 
-0.4750 
0.1530 




-0.9968 
-0.1150 
0.3880 
0.0850 



0.0077 
0.1430 




0.0802 0.0415 
-0.0318 
-0.4650 

1.0000 



(34) 



C 



10 
1 



For this example, n — 4, m = p = 2, so n = mp when q = 0, which is the 
dimension zero case. We can find some feedback laws with PHCpack. When 
the chosen poles are (—0.234, —1 + 3.2i, —1 — 3.2i, —3.0), we found two complex 
feedback laws. The total time spent is 140 milliseconds and the realization 
part just needs trivial time. By computing the eigenvalues of the closed-loop 
system with the MATLAB script, we find the difference between the computed 
eigenvalues and the given eigenvalues is just 10 - and the condition number 
is around 10. Although the result seems quite nice, we still want to find some 
real feedback laws. We could choose q = 1 to find some dynamic feedback laws, 
therefore n + q = 5 < mp + q(m + p) = 8, then we have three degrees of freedom. 
With the calculation of PHCpack and C interface, we get 8 solutions, including 
2 or 4, or 6 real solutions, depending on different additional eigenvalues. For 
a choice of (-0.234,-1 + 3.2i, — 1 - 3.2i, -3.0, -7.0), in which the first four 
are the same as in the static case, and three randomly generated poles are 
(—0.944,-0.995,-0.904), the relative difference of the computed eigenvalues 
vs. the given eigenvalues is bounded by 10 -10 . The total CPU time is 2 seconds 
and 680 milliseconds. For this choice of poles, we found four real feedback laws. 
For most of the 8 solutions, the condition numbers of the closed-loop system are 
less than 10 2 , and few of the solutions have the condition number 10 4 , calculated 
with formula (|31() . 



7 Conclusions 

In this paper we showed the practical feasibility of computing dynamic feedback 
laws using numerical homotopy algorithms, and we applied our software to 
examples from the control literature. Our publicly available implementation of 
the numerical greatest common divisor and numerical Smith normal form may 
be of independent interest. 
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